Propagation of hydrodynamic interactions between particles in a compressible fluid 
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Hydrodynamic interactions are transmitted by viscous diffusion and sound propagation: the tem- 
poral evolution of hydrodynamic interactions by both mechanisms is studied by direct numerical 
simulation in this paper. The hydrodynamic interactions for a system of two particles in a fluid 
are estimated by the velocity correlation of the particles. In an incompressible fluid, hydrodynamic 
interactions propagate instantaneously at the infinite speed of sound, followed by the temporal 
evolution of viscous diffusion. On the other hand, in a compressible fluid, sound propagates at a 
finite speed, which affects the temporal evolution of the hydrodynamic interactions by the order of 
magnitude relation between the time scales of viscous diffusion and sound propagation. The hydro- 
dynamic interactions are characterized by introducing the ratio of these time scales as an interactive 
compressibility factor. 



I. INTRODUCTION 

In particle dispersions, the motion of each particle in 
a fluid solvent affects the motion of the other particles. 
Such dynamical interactions are called hydrodynamic in- 
teractions; these interactions produce complex dynamical 
behavior for dispersions, which can be observed in parti- 
cle aggregation and sedimentation phenomena. 

Here, hydrodynamic interactions correspond to mo- 
mentum exchange among particles through the ambient 
fluid. Hydrodynamic interactions are transmitted by two 
mechanisms: viscous diffusion and sound propagation. 
These two mechanisms occur at different time scales. The 
time scale for viscous diffusion over a distance equal to 
the particle size is t u = a 2 /v, while the time scale for 
sound propagation is t c = a/c, where a is the particle 
radius, v is the kinematic viscosity, and c is the speed of 
sound in the fluid. To study dynamical effects at the par- 
ticle size scale, the relative significance of the sound prop- 
agation mechanism in hydrodynamic interactions can be 
assessed from the ratio of the two time scales described 
above: 



e - — - — 



(1) 



This dimensionless quantity is called the compressibility 
factor. Because sound propagation is much faster than 
viscous diffusion, the compressibility factor is generally 
quite small; for example, in a dispersion of water and 
particles of radius a = 100 nm, the compressibility fac- 
tor is estimated to be £ w 7 x 10~ 3 . In this case, the 
assumption of incompressibility is fully justified. Thus, 
fluids are assumed to be incompressible in many theo- 
retical studies of hydrodynamic interactions and viscous 
diffusion is only considered in temporal evolution mecha- 
nisms for hydrodynamic interactions 0, Q . On the other 
hand, the speed of sound in a liquid is around 10 3 m/s 
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irrespective of the liquid considered, producing a large 
compressibility factor for a dispersion of a highly viscous 
fluid solvent; for example, e sa 0.7 for olive oil and e ~ 10 
for corn syrup. 

In recent years, the temporal evolution of hydrody- 
namic interactions between two particles has been di- 
rectly observed 0-0]. In these experimental studies, par- 
ticles were trapped by optical tweezers and the correla- 
tions between the positional fluctuations of particles were 
measured. The authors reported the temporal evolu- 
tion of hydrodynamic interactions in the viscous diffusion 
regime, which coincided with analytical predictions that 
assumed the fluid was incompressible. Evidence for hy- 
drodynamic interactions by sound propagation was also 
observed [H 0, H| , but could not be captured due to the 
extremely short time scale of sound propagation. 

In the present study, we investigate the temporal evo- 
lution of hydrodynamic interactions by modeling sound 
propagation using direct numerical simulation. Within 
this approach, the hydrodynamic interactions are di- 
rectly computed directly by simultaneously solving for 
the motions of the fluid and the particles with appro- 
priate boundary conditions. We use the smoothed pro- 
file method (SPM) d QjJ, which can be applied to a 
compressible fluid as well as an incompressible fluid (llj . 
Sound propagation can only be captured by accounting 
for fluid compressibility. 

We consider a system of two particles in a fluid and 
investigate the correlated motion of the particles. In par- 
ticular, we estimate the velocity relaxation of the particle 
in response to the exertion of an impulsive force on the 
other particle. This velocity cross-relaxation function is 
equivalent to the velocity cross-correlation function from 
the fluctuation-dissipation theorem. The numerical re- 
sults for the velocity cross-relaxation function may be in- 
terpreted in terms of the temporal evolution of the flow 
field around the particle. We first consider an incom- 
pressible fluid to identify the separate characteristics of 
sound propagation and viscous diffusion in the hydro- 
dynamic interactions. We then consider a compressible 
fluid to consider the effect of the order of magnitude rela- 
tion of the sound propagation and viscous diffusion time 



scales on the temporal evolution of the hydrodynamic 
interactions. In addition, we examine the validity of an- 
alytical solutions within the Oseen approximation, which 
are often compared with experimental results. 

II. MODEL 

A. Basic equations 

We obtain solutions to describe the dynamics of parti- 
cles in a Newtonian fluid. The motion of the particles is 
governed by Newton's and Euler's equations of motion, 
which can be written for the i-th particle as 

Mi±Vi = Ff +Ff + F l E , ^Ri = Vi, (2) 



Particle 1 Particle 2 




FIG. 1. Geometry of the present model system. Two spherical 
particles of equal radius a are situated in a fluid with a center- 
to-center distance R = \R\. This system is axisymmetric 
about the R = R2 — R1 direction, such that two directions are 
of interest: parallel and perpendicular to the center-to-center 
axis of the particles. In this figure, the y and z directions are 
degenerate perpendicular directions. 



dt 



0,; = N, 



H 



(3) 



where Ri, V{, and i~2^ are the position, the translational 
velocity, and the rotational velocity of the i-th particle, 
respectively. The particle has a mass Mi and a moment 
of inertia J,. The fluid exerts a hydrodynamic force Ff 1 
and a torque Nf 1 on the particle, while a force Ff is ex- 
erted through direct interactions among the particles. A 
force F E and torque Nf are externally applied. The hy- 
drodynamic force and torque are evaluated from the mo- 
mentum conservation between the particle and the fluid. 

The fluid dynamics are governed by the following hy- 
drodynamic equations: 



| + v.( H = o, 

^ + V • (pvv) = V • er 
ot 



(4) 



(5) 



where p(r, t) and v(r,t) are the mass density and velocity 
fields, respectively, of the fluid. The stress tensor is given 
by 



-pi + rj[Vv + {Vv) T ] 



Vr 



n (v ■ v)i, 



(0) 



where p(r, t) is the pressure, r\ is the shear viscosity, and 
r\ v is the bulk viscosity. A body force f R (r, t) is added to 
satisfy particle rigidity. We also assume a barotropic fluid 
described by p = p(p) with a constant speed of sound c 
such that 



dp 



c 2 . 



(7) 



Equations dU)-© are closed for the variables p, v, and p, 
without consideration of energy conservation. 

We use SPM to perform direct numerical simulations. 
In this method, the boundaries between particle and fluid 
are modeled by a continuous interface. For this purpose, 



a smoothed profile function <f>(r,t) £ [0, 1] is introduced 
to distinguish between the particle and fluid domains: 
i.e., <j) = 1 in the particle domain and = in the 
fluid domain. These two domains are smoothly connected 
through a thin interfacial region of thickness £. The body 
force for the particle rigidity f R is expressed as p<j>f p . 
The detailed mathematical expressions of <f> and <f>f p have 
been previously given Q. 



B. Linear formulation 

The hydrodynamic equations can be linearized for slow 
flows in a dispersion. Correspondingly, the equations of 
motion of the particles are also linear. Most approximate 
theories are based on this formulation. 

Here, for a Ap-particle system, we define the following 
column vectors to describe the equations concisely: 



U 



V Np 



H 



*N P 

Nf 1 



( Fi E \ 



jpE 



which are the 6A p -dimensional vectors. Neglecting direct 
particle interactions, the equations of particle motion, 
i.e., Eqs. ([2|) and (J3]) are summarized by 



M--U(t)=H(t)+E(t), 



(8) 



where the matrix M represents the direct sum of the 
mass tensor and the moment of inertia tensor. For spher- 
ical particles with equal radii and densities, the matrix 
M reduces to 



M = 



MI 3Np 








3N„ 



(9) 
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where M and a are the mass and the radius, respectively, 
of each particle. With Fourier transforming Eq. ©, the 
corresponding equation for the Fourier components with 
the time factor e~ lut is obtained by 



iuM ■ U{uj) = H{u) + E(lu). 



(10) 



The following relation holds true for the linearized hy- 
drodynamic equations Q: 



U(u) = • H(oj), 



(11) 



where A is the 6N P x 6N P mobility matrix. The mobil- 
ity matrix depends on the configuration of all the parti- 
cles {Ri}, so we further assume, for simplicity, that the 
particles do not move significantly over the time scale 
considered, namely, the configuration of particles {Ri} 
is independent of time t. Substitution of Eq. into 
Eq. (|10p yields an explicit expression for the translational 
and rotational particle velocities: 

U(u>) = [I - iujp,(u) • M]" 1 • AM ■ E(cj). (12) 

Therefore, for an externally applied force and torque, the 
dynamics of particles in a fluid can be described by the 
mobility matrix. 

The components of the mobility matrix are given by 



A(w) 



A U M A tr H 
A r » A rr M 



(13) 



where A Q/3 is the 3N p x 3N p matrix, which is composed 
of the mobility tensors Ay for two particles i and j. 
The mobility tensors describe the mutual coupling of 
the translational and rotational motion between parti- 
cles and the superscript represents the mode of the mo- 
tion: 't' for translation and 'r' for rotation. The mobility 
matrix is symmetric as a consequence of the Lorentz re- 
ciprocal theorem, such that the following relations are 
satisfied 0: 



(A 



tt\T 



rr\T 



(A rr ) 



(A 



rt\T 



(14) 



Now, we consider a system of two identical spherical 
particles in a fluid. The configuration of the particles is 
described in Fig. [TJ Numbers are assigned to the parti- 
cles: 1 for the particle on the left and 2 for the particle on 
the right. We define a vector R — R2—R1 to describe the 
geometry of this system. The particle center-to-center 
distance is denoted by R = \R\. The unit vector along 
the line of centers is described by R — R/R. Due to 
axisymmetry about the R axis, each mobility tensor is 
described by at most two scalar functions [ill 13 : 



A*f (R,u)(I-RR), 

(15a) 



u tr± 



{R,lu)Rx I. 



(15c) 



The superscripts || and _L denote the directions of motion 
parallel and perpendicular to the symmetry axis, respec- 
tively. In the direction parallel to the symmetry axis, the 
translational and rotational motions are decoupled as in- 
dicated by Eq. (|15cl) . Interchanging the particle numbers 
1 and 2 corresponds to inverting the direction of R, which 
causes sign inversion for the translation-rotation cross- 
mobility tensor Ay without changing sign of the other 
cross- mobility tensors Ay an d Ay • The mobility tensors 
Eqs. (I15p within the Oseen approximation are described 
in the Appendix. 



C. Temporal evolution of the flow field 

We first consider a particle dispersion with an incom- 
pressible fluid solvent. Incompressibility corresponds to 
an infinite speed of sound in the fluid, such that the 
solenoidal condition is imposed on the fluid velocity field 
from Eqs. (gj) and © as 



0. 



(16) 



The velocity field can be decomposed into two contribu- 
tions as 



w — Vy. 



(17) 



The vector field w(r,t) is the juxtaposition of the ve- 
locity fields in the fluid and particle domains, so that 
the fluid-particle impermeability boundary condition, or 
the continuity of the normal velocity on the boundary, 
is not satisfied. An irrotational flow field, described by 
the scalar potential <p(r,t), is added in order that the 
velocity field v satisfies the boundary condition. Due to 
the solenoidal condition on the velocity field v, the scalar 
potential ip obeys the Poisson equation: 



V 2 (yS = V • W. 



(18) 



On the right-hand side, V • w is zero except at the fluid- 
particle boundary, at which singularities exist. From 
Eq. (jT5J), the scalar potential is given by [1 5j j 



<p(r,t) 



1 

47T 



dr 



,{r- 



r — r 



l\3 



w{r',t). (19) 



Now we consider an isolated single particle, due to ax- 
isymmetry, the scalar potential and the corresponding 
velocity field can be simply described as 



Q(t) 



Vcp(r,t) 



Q(t) ■ (3rr - I) 



(20) 



(21) 



A«(-H.w) 



f$ l (R,Lo)RR- 



Alf (i^X-T - RR), 



(15b) 



where r — r jr is a unit vector in the radial direction. The 
velocity field of Eq. (j2"Tj) represents a doublet flow, which 
corresponds to the electric field generated by an electric 
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dipole. The vector Q(t) is parallel to the particle velocity 
and depends on time in a reflection of the particle motion 
relative to the fluid. The analytical form of Q(t) has been 
derived previously [16|, and the strength Q(t) = \Q(t)\ is 
time-independent for a neutrally buoyant particle. The 
doublet flow is generated instantaneously to satisfy the 
solenoidal condition, thus, the doublet flow is interpreted 
to expand at the infinite speed of sound. 

The vector field w represents the shear flow due to 
viscous diffusion. For sufficiently slow fluid flow, the hy- 
drodynamic equations can be linearized as follows: 



dv 



7 v + f 



(22) 



Under the incompressibility condition, the pressure gra- 
dient imposes fluid-particle impermeability on the time- 
derivative of the velocity field, so that the pressure is 
related to the scalar potential as 



dip 



(23) 



Therefore, the temporal evolution of the vector field w is 
given by 



dw 2 __ 

at 



w) 



(24) 



Modifying the diffusion term on the right-hand side to 
be solenoidal, results in the doublet flow being gradually 
offset by the diffusion of the vector field w . 

In short, the exertion of an impulsive force on a parti- 
cle in an incompressible fluid instantaneously generates 
a doublet flow. The shear flow then spreads out dif- 
fusely, cancelling out the doublet flow. This picture will 
be validated in one of the following numerical simulation 
results. 

For a compressible fluid, a second contribution is added 
to Eq. (HZl): 



v = w — V-0 — V<^, 



(25) 



where the scalar potential i/j(r,t) corresponds to the bulk 
compression flow. Here, we assume that the dynamics of 
fluid can be described by the linearized hydrodynamic 
equations: 



dp 

at 



-p c V • v, 



(26) 



dv { 1 \ 

Po ^=7]V 2 v+ i-^V + Vv) W-v-Vp + f R (27) 

The time evolution equation of the total scalar potential 
tp' = tp + if) is derived from the equations above as 



vv 



a , 1 0V 



c 2 dt 2 p c 2 dt 



= I 



Poc 



2 V 2 -)(V-™), (28) 



where iji — (4/3)t] + t] v denotes the longitudinal viscosity. 
This is a damped wave equation, where the source is 
given on the right-hand side of the equation. Therefore, 
a potential flow propagates at a finite speed of sound 
in a compressible fluid. In the long-time limit, Eq. (|28|) 
reduces to the Poisson equation for an incompressible 
fluid, as given by Eq. (TT5|) . and the potential flow reduces 
to an instantaneous doublet flow. 



III. NUMERICAL RESULTS 

Numerical simulations are performed for a three- 
dimensional box with periodic boundary conditions. The 
space is divided into meshes of length A, which is a unit 
length. The units of the other physical quantities are de- 
fined by combining 77 = I and po = I with A, where po is 
the fluid mass density at equilibrium. The system size is 
L x x L y x L z = 256 x 256 x 256. The other parameters 
are set to a = 4, £ = 2, p p = 1, rj v = 0, and h = 0.05, 
where p p is the particle mass density and h is the time 
increment of a single simulation step. 

We consider a system in which two identical spherical 
particles in a fluid, whose geometry is described in Fig[T] 
Both particles have the same density and radius, namely 
the same mass and moment of inertia. We investigate 
the time-dependence of the velocity for particle I fol- 
lowing the exertion of an impulsive force at the center of 
particle 2, with changing the center-to-center distance be- 
tween particles given by R* = R/a. This cross- relaxation 
function is a manifestation of the temporal evolution of 
the hydrodynamic interactions between particles I and 
2. The impulsive force is assumed to be sufficiently small 
such that the Reynolds and Mach numbers of the flow are 
sufficiently low. We set the impulsive force to produce 
an initial particle Reynolds number of Re p = I0~ 3 . In 
addition, because the particle displacement before stop- 
ping, scaled by the particle radius a, is comparable to 
{p P l po)Ke p , the displacement of particles is negligible 
in the present simulations. Therefore, because particle 
collisions do not occur, the direct interactions between 
particles, including overlap repulsion forces, are not con- 
sidered. 

The cross-relaxation tensor is introduced as the nor- 
malized change in velocity of particle 1: 



V 1 (R 1 ,f < R 2 ) = ^-7i2(R,t), 



(29) 



ji 2 (R,t) ^ 1 [ 2 {R,t)RR + tf 2 (R,t)(I - RR), (30) 

where P% is the impulsive force exerted on particle 2 at 
t = 0. Due to the previously mentioned axisymmctry, 
the cross-relaxation tensor is characterized by only two 
directions of motion, which are parallel and perpendicu- 
lar to the center-to-center vector R. Then, the motions 
in both directions are decoupled. From the fluctuation- 
dissipation theorem, the cross-relaxation tensor is equiv- 
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alent to the velocity cross-correlation function of the fluc- 
tuating system: 



712 (R,t) = 



M 



;vi(o,o)v 2 (jm)), 



(31) 



where ks is the Boltzmann constant and T is the ther- 
modynamic temperature. We examine both the parallel 
correlation jf 2 an( l the perpendicular correlation by 
adjusting the direction of the impulsive force P%. As 
predicted by the formulation in the preceding section, 
for motion perpendicular to the axis, the translation- 
rotation coupling in the particle motion is also observed. 
However, because the influence of sound propagation on 
the rotational motion is expected to be small, we will fo- 
cus entirely on the translation-translation coupling. We 
also compare the simulated results with approximate so- 
lutions, where the analytical solution for an isolated sin- 
gle particle for the self-mobilities and the Oseen approx- 
imation for the cross-mobilities are applied (see the Ap- 
pendix) . 



A. Incompressible fluid 

First, we investigate the temporal evolution of hydro- 
dynamic interactions in an incompressible fluid. The dy- 
namics of the incompressible fluid is governed by the hy- 
drodynamic equations composed of Eqs. |5]), ©, and 
(fT6|) . The center-to-center distance between particles 
takes values of R* — 3, 4, and 7. The simulation re- 
sults for the velocity cross-relaxation functions are shown 
in Fig. [21 While the parallel correlation ~i\ 2 {t) is pos- 
itive at all times, The perpendicular correlation 7f^(i) 
is initially negative and subsequently becomes positive. 
Comparable results have been reported in experimental 
studies 0-IH], where the difference between the parallel 
and perpendicular correlations has been attributed to the 
flow field around the particle [J, 0] ■ We discuss the time- 
dependence of the velocity cross-relaxation functions in 
more detail below. 

The temporal evolution of the flow field has been ex- 
plained in the previous section. Here, the velocity fields 
generated by the impulsive force exerted on a single par- 
ticle are shown in Fig. [3] These simulation results are ob- 
tained for a single particle system with size L x xL y xL z = 
128 x 128 x 128. At early times t/r v = 2.5 x 1CT 2 , the dou- 
blet flow is dominant. The doublet flow is characterized 
by loop streamlines flared in the direction perpendicu- 
lar to the particle motion, where backflow is observed as 
described by Eq. (|20[) . As discussed in the previous sec- 
tion, the doublet flow appears instantaneously, which is 
interpreted as the propagation of an infinite-speed sound 
wave. On the other hand, the shear flow due to viscous 
diffusion, whose strength is given by the vorticity Vxv, 
is only observed in the vicinity of the particle. Over a 
time t/T„ =1.2, the shear flow has diffused over a large 
range following Eq. (|2"4")l ; therefore, there is a spreading 
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FIG. 2. Velocity cross-relaxation function in an incompress- 
ible fluid in (a) parallel and (b) perpendicular directions rel- 
ative to the symmetry axis of the particles. The distance 
between the two particles are R* = 3, 4, and 7. Simulation 
results are described by bold solid lines (positive values) and 
dashed-dotted lines (negative values) . Solutions within the 
Oseen approximation are represented by thin solid lines (pos- 
itive values) and dashed-dotted lines (negative values). The 
broken line represents the long-time tail with the algebraic 
power law decay Bt~ 3 ^ 2 , as given by Eq. 



fluid region flowing in the same direction as the parti- 
cle motion, while the loop streamlines get away from the 
particle. 

The velocity cross-relaxation function is related to the 
temporal evolution of the flow field around particle 2. 
The doublet flow propagated by sound waves produces 
an instantaneous velocity correlation between the parti- 
cles. The doublet flow direction shown in Fig.[3]produces 
the positive parallel correlation and the negative perpen- 
dicular correlation. Little change in the cross-relaxation 
functions occurs in the early stages, reflecting the time- 
independence of the strength of the doublet flow due to 
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(a) t/x v = 2.5 xlO" 2 (b) t/x v = 1.2 




FIG. 3. Temporal evolution of the velocity field around an isolated single particle due to an impulsive force exerted on the 
particle in an incompressible fluid. The cross-sections parallel to the impulsive force direction including the particle center are 
shown. The impulsive force is exerted at time t — and the simulation results are given at (a) t/r u = 2.5 x 10~ 2 and (b) 
t/Tu = 1-2. The direction of the impulsive force is to the right in the pictures and the particle is represented by a black circle. 
The velocity fields are plotted only for |i>|(a/^Re p ) > 10~ 4 . The vorticity of the velocity field V x v is described by a color 
scale (gray scale), going from negative (darker) to positive (lighter) vorticity. The vorticity is normalized by T„/Re p . 



(a) t/x v = 6.25 xlO" 1 (parallel) 
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(c) t/x v = 12.5 (parallel) 




(b) t/x v = 6.25 xlO" (perpendicular) 
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(d) t/x v = 12.5 (perpendicular) 
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FIG. 4. Temporal evolution of the velocity field around the particles. The cross-sections including the center-to-center axis 
between particles are shown. The center-to-center distance between two particles is R* — 7. The impulsive force is exerted at 
time t — on the particle on the right and the simulation results are presented at (a,b) t/r u = 6.25 x 10~ and (c,d) t/r v = 12.5. 
The direction of the impulsive force is (a,c) to the right and (b,d) upwards. The particles are represented by black circles. The 
velocity fields are shown only for |v|(a/i^Re p ) > 1CP 4 . The vorticity of the velocity field V x v is described by a color scale 
(gray scale), going from negative (darker) to positive (lighter) vorticity. The vorticity is normalized by T„/Re p . 
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FIG. 5. Velocity cross-relaxation function in a compressible fluid in (a, b) parallel and (c, d) perpendicular directions to the 
symmetry axis. The center-to-center distance between two particles is R* = 3. The compressibility factor takes values of 
e = 0.1, 0.6, and 1.5. Simulation results are described by bold solid lines (positive values) and dashed-dotted lines (negative 
values). Solutions within the Oseen approximation are represented by thin solid lines (positive values) and dashed-dotted lines 
(negative values). 



the neutrally buoyant particle [16j . 

The reduction in the parallel correlation and the sign 
inversion of the perpendicular correlation occur at about 
the same time, which depends on the inter-particle dis- 
tance R*. Because shear flow by viscous diffusion can 
cause such changes in the cross-relaxation functions, we 
estimate the time scale on which the shear flow gener- 
ated by the particle 2 arrives at particle 1 by the viscous 
diffusion time scale over the length L: t* = L 2 /v. In this 
case, the characteristic length is the distance between the 
particle surfaces, L = R — 2a; thus, for R* — 3, 4, and 7, 
the viscous diffusion time scales are t* fr v = (R* — 2) 2 = 
1, 4, and 25, respectively. The time for the parallel cor- 
relation to begin decreasing and for the perpendicular 
correlation to invert the sign roughly corresponds to r* 
in each case. 



In an incompressible fluid, hydrodynamic interactions 
are instantly propagated at the infinite speed of sound 
and subsequently transmitted by viscous diffusion in a 
time scale t* . The temporal evolution of the flow fields 
around two particles separated by R* = 7 is shown in 
Fig. |U At early times in the flow, t/r v — 6.25 x 1CP 1 so 
that the shear flow does not arrive at particle 1, whose 
dynamics are governed by the doublet flow. Later, at 
t/r u = 12.5, the shear flow diffuses over time to arrive at 
particle 1 and the motion of particle 1 is governed by the 
shear flow. Here, the stresslet flow field is observed in the 
vicinity of particle 1. This flow is generated to maintain 
particle rigidity against deformation in a shear flow. In 
the perpendicular correlation, the rotational motion of 
particle 1 is also observed. 

In the final stages of the flow, the momentum of parti- 
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FIG. 6. Velocity cross-relaxation function in a compressible fluid in (a, b) parallel and (c, d) perpendicular directions to the 
symmetry axis. The center-to-center distance between two particles is R* = 7. The compressibility factor takes values of 
e = 0.1, 0.6, and 1.5. Simulation results are described by bold solid lines (positive values) and dashed-dotted lines (negative 
values). Solutions within the Oseen approximation are represented by thin solid lines (positive values) and dashed-dotted lines 
(negative values). 



cle 2 has completely diffused away, so that the particles 
and the fluid move collectively. These dynamics are man- 
ifested in a long-time tail with a power law decay i~ 3 / 2 
of the cross-relaxation functions shown in Fig. [5J which 
are given by [ItJ 

7^(0 = -^^(^) 3/2 as t->oo. (32) 

Discrepancies between the simulation results and the 
Oseen approximation are observed, especially at t/r v < 
3, as shown in Fig. [2j The approximate cross-relaxation 
functions clearly change with time in the early stages and 
the time at which the effect of the shear flow becomes 
apparent is later than in the simulation results. Parti- 
cles are assumed to be points in the Oseen approxima- 
tion; therefore, this approximation implies that the par- 
ticle separation is much larger than the particle radius as 



R* = JJ/b> 1 and that the time scale of observation is 
much longer than that of viscous diffusion over the length 
of the particle radius as tjr v = tv/a 2 ^> 1. Neglecting 
the particle size results in a larger characteristic length 
L = R and a longer diffusion time scale t*/t„ = R* 2 . 
Consequently, the Oseen approximation can only be ac- 
curately applied at long times and large particle sepa- 
rations. The simulation results shown in Fig. [2] confirm 
that the validity of the Oseen approximation increases 
with increasing inter-particle distance R*. In previous 
experimental studies, measurements were conducted over 
the range for which the Oseen approximation is valid, 
i.e., R* > 5 and measuring frequencies corresponding to 
tjr v > 30: consequently, the experimental results showed 
good agreement with the Oseen approximation 0, 
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(a) 8 = 0.1 (parallel) (b) 8 = 0.1 (perpendicular) 




(c) 8 = 1.5 (parallel) (d) 8 = 1.5 (perpendicular) 




FIG. 7. Velocity field around particles in a compressible fluid. Simulation results at t/r v — 6.25 x 10" 1 for compressibility 
factors with values (a,b) e = 0.1 and (c,d) e = 1.5. The velocity fields are shown only for \v\(a/v~Re p ) > 10 -4 . The divergence 
of the velocity field V v is described by a color scale (gray scale), going from negative (darker) to positive (lighter) divergence. 
The value of divergence is normalized by r^/Ke p . 



B. Compressible fluid 

Next, we consider a compressible fluid in which 
sound propagates at a finite speed. The velocity cross- 
relaxation functions for inter-particle distances R* = 3 
and 7 are shown in Figs. \S\ and [6j respectively. Due to 
the periodic boundary conditions imposed on the simula- 
tion box, a sound pulse generated by the particle motion 
can affect the simulation results after the sound pulse ar- 
rives at the end of the simulation box. The oscillatory 
structure of the cross-relaxation functions, especially at 
e = 0.1, is one of the striking artifacts of the periodic 
boundary conditions. 

The velocity cross-relaxation functions are zero in the 
early stages and then suddenly change to non-zero values. 
Because this time lag increases with the inter-particle 
separation, this is likely to be due to the sound propaga- 
tion from particle 2 to particle 1 . As with the diffusion of 
shear flow, we estimate the time scale for sound propaga- 
tion over a characteristic length L — R — 2a as r* = L/c. 
For R* = 3 and 7, the sound propagation time scales 
are t*/t v = (R* — 2)e = e and 5e, respectively. The 



results shown in Figs. [5] and [5] confirm the correspon- 
dence between the time scale of sound propagation and 
the peak of the cross-relaxation functions. The peak in 
the cross-relaxation function broadens as sound propaga- 
tion is dissipated by the longitudinal viscosity according 
to Eq. (EH]). 

When the fluid compressibility is as small as e — 0.1, 
the velocity cross-relaxation functions superimpose onto 
those for the incompressible fluid after a time r*. This 
behavior only describes hydrodynamic interactions by 
sound propagation before viscous diffusion effects come 
into play. On the other hand, for a large compressibility 
e = 1.5, the behavior of the cross-relaxation functions 
are considerably different, except for the hydrodynamic 
long-time tail: in particular, there is no negative perpen- 
dicular correlation. For a fluid with medium compress- 
ibility such as e = 0.6, the negative perpendicular cor- 
relation is observed only for inter-particle separations as 
large as R* — 7. As sound propagation and viscous dif- 
fusion generate negative and positive perpendicular cor- 
relations, respectively, the balance of the two time scales 
t* and r* is expected to characterize the behavior of the 



10 



(a) divergence (parallel) 



(b) divergence (perpendicular) 




11 llllllll 

(c) vorticity (parallel) 




(d) vorticity (perpendicular) 




FIG. 8. Velocity field around particles in a compressible fluid. Simulation results at t/r v = 3.13 for compressibility factor 
with value s — 1.5. The velocity fields are shown only for \v\(a/ uTLe p ) > 10~ 4 . The (a,b) divergence and (c,d) vorticity of 
the velocity field are described by a color scale (gray scale), going from negative (darker) to positive (lighter) divergence. The 
value of divergence is normalized by T„/Re p . 



cross-relaxation function. Therefore, we define the inter- 
active compressibility factor using a characteristic length 
L = R — 2a as 



As the inter-particle distance increases, the two time 
scales separate to reduce the interactive compressibility 
factor. In the present simulations, the interactive com- 
pressibility factors for R* — 3 and 7 are e* — e and 0.2e, 
respectively. When the interactive compressibility factor 
is small, hydrodynamic interactions propagating at the 
speed of sound arrive first at the other particle in about 
a time r* , followed by those propagating by viscous dif- 
fusion arriving in about a time r*. The temporal evo- 
lution is qualitatively the same as in an incompressible 
fluid, except for the time lag due to sound propagation 
at a finite speed. This situation applies to the results 
where a negative perpendicular correlation is observed. 
On the other hand, if the interactive compressibility fac- 
tor is sufficiently large, the order of arrival of sound prop- 
agation and viscous diffusion should be reversed. In this 
situation, the perpendicular correlation is positive from 



the start due to viscous diffusion and the effect of sound 
propagation on the hydrodynamic interactions produces 
a local minimum at most in the cross-correlation func- 
tion. In the present studies, such results are observed 
when at least e* > 0.3 is satisfied. 

The velocity fields around particles separated by R* = 
7 at i/r„ = 6.25 x 10 _1 are shown in Fig. [7] The expand- 
ing doublet flow with a strength given by the divergence 
of the velocity field V • v are observed. For a small com- 
pressibility factor e = 0.1, the doublet flow has just ar- 
rived at particle 1. This situation corresponds to a peak 
in the cross- relaxation functions shown in Fig.[5ta,c). On 
the other hand, for a large compressibility factor e = 1.5, 
the doublet flow has not yet arrived at particle 1 and the 
cross-relaxation functions are correspondingly zero. At a 
later time t/r u — 3.13, as shown in Fig. [5J the shear flow 
(vorticity) arrives at particle 1, instead of the doublet 
flow (divergence) . These dynamics are manifested in the 
positive cross-relaxation functions of the perpendicular 
correlation shown in Fig. |Hl[b,d). 

In a compressible fluid, because there is a finite sound 
propagation time scale as well as a viscous diffusion time 
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scale, the time condition for the validity of the Oseen 
approximation is given by t 3> max(r 1/ ,r c ). As with the 
diffusion of shear flow, sound propagation from particle 2 
to particle 1 is delayed within the Oseen approximation 
as shown in Figs. [5] and [6] This is because the sound 
propagation time scale is increased within this approx- 
imation so that t*/t„ = R*e. Thus, the discrepancies 
with the simulation results can be fairly well-resolved by 
increasing the particle separation. 



ACKNOWLEDGEMENTS 

This work was supported by KAKENHI 23244087 and 
the JSPS Core-to-Core Program "International research 
network for non-equilibrium dynamics of soft matter." 

APPENDIX: MOBILITY MATRIX IN THE 
OSEEN APPROXIMATION 



IV. CONCLUSION 

In the present study, we investigated the temporal evo- 
lution of hydrodynamic interactions for a system of two 
particles, using SPM to perform direct numerical simu- 
lations. Fluid compressibility was considered in examin- 
ing temporal evolution by sound propagation. Hydrody- 
namic interactions were estimated by the velocity cross- 
relaxation functions, which are equivalent to the velocity 
cross-correlation functions in a fluctuating system. 

In an incompressible fluid, hydrodynamic interactions 
were observed to propagate instantaneously at the infi- 
nite speed of sound, along with subsequent temporal evo- 
lution by viscous diffusion. Theoretical analysis showed 
that a doublet flow and shear flow are generated by sound 
propagation and viscous diffusion, respectively. Because 
the cross-relaxation function reflects the characteristics of 
each flow field, sound propagation and viscous diffusion 
are associated with negative and positive perpendicular 
correlations between the partice velocities, respectively. 
Therefore, the time at which the behavior of the cross- 
relaxation function changes is related to the time scale of 
viscous diffusion over the particle separation r* , which is 
only the time scale characterizing the temporal evolution 
of hydrodynamic interactions in an incompressible fluid. 

In a compressible fluid, sound propagates between par- 
ticles in a finite time, which scales as r* . The effect of the 
order of magnitude relation between the two time scales 
t* and t* on the cross-relaxation function was observed; 
an interactive compressibility factor e* was defined as 
the ratio of the two time scales given by Eq. (|33[) . In 
our simulation results, the reversal of the order of ar- 
rival of sound propagation and viscous diffusion at the 
other particle was observed at e* > 0.3, for which there 
was no negative perpendicular correlation. In this case, 
the hydrodynamic interactions were largely governed by 
viscous diffusion. 

The temporal evolution of hydrodynamic interactions 
does not qualitatively change from that for an incom- 
pressible fluid, as long as the interactive compressibility 
factor is small. Only when the interactive compressibility 
factor is sufficiently large can differences in the temporal 
evolution of hydrodynamic interactions be expected. As 
such a situation can be realized in a highly viscous fluid, 
experimental validation of the results presented here is 
desirable. 



Hydrodynamic interactions among particles in slow 
flows are described by the mobility matrix discussed in 
§2.2. However, there is no known closed-form solution 
even for a two-particle system. In this section, we in- 
troduce an approximate analytical form for the mobility 
tensors. 

For slow flows, the hydrodynamic equations can be 
linearized as given by Eqs. ([2^)1 and ([2"T]) . and the cor- 
responding equations for the Fourier components, with 
the time factor e~ luJt , are obtained as 



Poc 2 V • v u = 0, 



(Al) 



- iu)v u = rjV 2 ^ + i-r) + r) v j VV • v u - Vp u + f*. 

(A2) 

The fluid velocity field generated by the body force f R 
is expressed as 



v u (r) = / dr'G(r-r» • /*(r') 



(A3) 



where the Green's function is given by [18j, 1 1 9j j 

pijlT _ p -ar 



1 / e~ a 
4ttt] \ r 



with 



and 



a = (-iujp /i]) 1/2 , h = uj/c, 



c = c 



iuj ( 4 
p c z V 3 



1/2 



(A4) 
(A5) 

(A6) 



For an isolated single spherical particle in a fluid, with 
a body force constraint f R to satisfy the condition of 
particle rigi dity , the self-mobilities can be analytically 
derived [ia|M|2llas follows: 



t 9 2x 2 (l - iy) - (1 + x)y 2 - x 2 y 2 



2x 2 (1 + x)(9 - 9iy - 2y 2 ) + x 2 (l - iy) ' 

(A7) 



„ rr , . r 3(1 + x) 



] 3 + 3x + x 2 ' 



(A8) 
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where x = aa, y = /ia, /1q = (67r?7a) , and /itg = 
(Swrja 3 )^ 1 . When the speed of sound is assumed to be 
infinite, or = 0, the solution for an incompressible fluid 
is obtained. The translation-rotation coupling does not 
appear in the self-mobility as Aii = 0- 

To calculate the cross-mobilities, we introduce the Os- 
een approximation in which the particles are regarded as 
points. In the Oseen approximation, a particle exerts a 
Stokeslet (point force); therefore, according to Eq. (|A3|) . 
the translation-translation cross-mobility is the Green's 
function itself: 

P? 2 (R,l>) = G(R,u,). (A9) 

The components of the parallel and perpendicular direc- 
tions are explicitly given by 

~(2 + 2aR)e- aR ], (AlOa) 

^(R,^ = 0* + i^R)e^ R 

+(l + aR + a 2 R 2 )e- aR }. (AlOb) 

The Stokeslet also generates a rotational motion with 
an angular velocity of (1/2)V x v u ; therefore, the 



translation-rotation cross-mobility is given by 

P*(R,cj) = -^V x G(R,uj), (All) 
which has a perpendicular component only: 

P&-(R,u) = ^(l + a R)e- aR . (A12) 

There are no effects from the fluid compressibility. This 
result is also derived by considering the velocity field 
generated by a rotlet (point torque) [22J. The rotation- 
rotation cross-mobility is obtained from the angular ve- 
locity generated by a rotlet: 

AS(^w) = \ v x v x G(R,u). (A13) 
Then, each component is given by 

AS" (R, w) = Mo|J(l + ocR)e- aR , (Al4a) 

fj&iRiu) = -M r o^s (1 + «i? + a 2 R 2 )e- aR . 

(A14b) 

Using the analytical solutions for an isolated single 
particle for the self-mobilities and the Oseen approxima- 
tion for the cross-mobilities, the velocity cross-relaxation 
functions can be approximated by Eq. Q12[), 
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